#include <stdlib.h>
#include <stdio.h>
#include "nbody.h"
#include <math.h>
#include <string.h>
#include <pthread.h>
#include <time.h>
#define MAX 255

int numCuerpos, numPasos, anchuraTablero, alturaTablero, numProcesos;
double g, intervalo, minDist;
double *posicion[2];
double *velocidad[2];
double *fuerza[2];
double *masa;
pthread_mutex_t mp;


int main(int argc, char **argv){
    int i=0;
    Hilos *hilo;
    
    numProcesos=(atoi(argv[3]));
    leerCabecera(argv[1]);
    
    for(i=0;i<numPasos;i++){
        crearHilosCalcFuerzas(hilo);
        crearHilosCalcPosicion(hilo);
    }
    printf("Guardando datos en el fichero: %s\n", argv[2]);
    escribirDatos(argv[2]);
    
    return 1;
}

    int crearHilosCalcFuerzas(Hilos *hilo) {
        int i=0;
        hilo=(Hilos *)malloc(numProcesos*sizeof(Hilos));
        for (;i<numProcesos;i++) {
            hilo[i].idHilo=i;
            
            if(pthread_create(&hilo[i].thread, NULL, calcularFuerzas,  hilo[i].idHilo))
                printf("Error en la creacion del hilo [%d]\n", hilo[i].idHilo);
        }
        for (i=0;i<numProcesos;i++)
            pthread_join(hilo[i].thread, NULL);    
  
        return 1;
    }
    
    int crearHilosCalcPosicion(Hilos *hilo) {
        int i=0;
        hilo=(Hilos *)malloc(numProcesos*sizeof(Hilos));
        for (;i<numProcesos;i++) {
            hilo[i].idHilo=i;
            
            if(pthread_create(&hilo[i].thread, NULL, calcularPosicion,  hilo[i].idHilo))
                  printf("Error en la creacion del hilo [%d]\n", hilo[i].idHilo);
     }
        for (i=0;i<numProcesos;i++) 
            pthread_join(hilo[i].thread, NULL);   
        
        return 1;
    }

int leerCabecera(char *fichero_entrada){
    FILE *file;
    int  cont;
    
    //Apertura del fichero y lectura de la cabecera
    
    file = fopen(fichero_entrada, "r");
    
    fscanf(file, "%d %lf %lf %d %lf %d %d", &numCuerpos, &g, &intervalo, &numPasos, &minDist, &anchuraTablero, &alturaTablero);
    
    inicializarDatos();
    
    printf("numcuerpos: %d, g: %lf, intervalo: %lf, numPasos: %d,mindist: %lf,anchura: %d,altura: %d\n", numCuerpos, g, intervalo, numPasos, minDist, anchuraTablero, alturaTablero);

    for(cont=0; cont<numCuerpos; cont++){
        fscanf(file, "%lf %lf %lf %lf %lf \n", &posicion[0][cont], &posicion[1][cont], &velocidad[0][cont], &velocidad[1][cont], &masa[cont]);
 }
  return 1;
}

int escribirDatos(char *fichero_salida){
    FILE *file;
    int cont;
    
    //Apertura del fichero y escritura de la cabecera
    file = fopen(fichero_salida, "wr");
    fprintf(file, "%d %lf %lf %d %lf %d %d\n", numCuerpos, g, intervalo, numPasos, minDist, anchuraTablero, alturaTablero);
    
    //Escritura de los datos
    for(cont=0; cont<numCuerpos; cont++){
        fprintf(file, "%lf %lf %lf %lf %lf\n", posicion[0][cont], posicion[1][cont], velocidad[0][cont], velocidad[1][cont], masa[cont]);
    }
    fclose(file);
    return 1;
}

void inicializarDatos(){
    posicion[0] = (double *)malloc(numCuerpos* sizeof (double));
    posicion[1] = (double *)malloc(numCuerpos* sizeof (double));
    fuerza[0] = (double *)malloc((numCuerpos)* sizeof (double));
    fuerza[1] = (double *)malloc((numCuerpos)* sizeof (double));
    velocidad[0] = (double *)malloc((numCuerpos)* sizeof (double));
    velocidad[1] = (double *)malloc((numCuerpos)* sizeof (double));
    masa = (double *)malloc(numCuerpos* sizeof (double));
}

void *calcularFuerzas(int proceso){
    double fax, fay, dist, f, alpha, aux;
    int i, j, planetaInicial, cuerposProceso;
    
   
    cuerposProceso=(numCuerpos/numProcesos);
    planetaInicial=cuerposProceso*proceso;
   
    for (i = 0; i < numCuerpos&&proceso==0; i++) {
        fuerza[0][i] = fuerza[1][i] = 0;
    }
  
    for (i = planetaInicial; i < planetaInicial+cuerposProceso; i++) {
        for (j = i+1; j < numCuerpos; j++) {     
            dist = sqrt((posicion[0][j] - posicion[0][i]) * (posicion[0][j] - posicion[0][i]) + (posicion[1][j] - posicion[1][i]) * (posicion[1][j] - posicion[1][i]));   
            if (dist > minDist ) {
                f = g * masa[i] * masa[j] / (dist * dist);
                
                if (f > MAX_F)
                    f = MAX_F;
                
                aux  = (posicion[1][j]-posicion[1][i]) / (posicion[0][j]-posicion[0][i]);
                
                alpha = atan(aux);
               
                fax   = f * mycos(alpha);
                fay   = f * mysin(alpha);
              
                pthread_mutex_lock(&mp);
                fuerza[0][i] += fax;
                fuerza[1][i] += fay;
                fuerza[0][j] -= fax;
                fuerza[1][j] -= fay;
                pthread_mutex_unlock(&mp);
                
             
            }
            
        }
    }

 pthread_exit((void*)0);
}

void *calcularPosicion(int proceso){
    double *ax, *ay;
    int i, planetaInicial, cuerposProceso;
    
    cuerposProceso=(numCuerpos/numProcesos);
    planetaInicial=cuerposProceso*proceso;

    ax = (double *)malloc(sizeof(double) * numCuerpos);
    ay = (double *)malloc(sizeof(double) * numCuerpos);

    for (i = planetaInicial; i < planetaInicial+cuerposProceso; i++) {
  
        //calculo de las aceleraciones
        ax[i] = fuerza[0][i] / masa[i];
        ay[i] = fuerza[1][i] / masa[i];
        //calculo de las velocidades
        velocidad[0][i] = velocidad[0][i] + ax[i] * intervalo;
        velocidad[1][i] = velocidad[1][i] + ay[i] * intervalo;
        //calculo de las posiciones
        posicion[0][i] = posicion[0][i] + velocidad[0][i] * intervalo;
        posicion[1][i] = posicion[1][i] + velocidad[1][i] * intervalo;
       
        if (posicion[0][i] <= 0) {
            posicion[0][i] = 2;
            velocidad[0][i] = - velocidad[0][i];
        }
        if (posicion[0][i] >= anchuraTablero ) {
            posicion[0][i] = anchuraTablero - 2;
            velocidad[0][i] = - velocidad[0][i];
        }
        if (posicion[1][i] <= 0) {
            posicion[1][i] = 2;
            velocidad[1][i] = - velocidad[1][i];
        }
        if (posicion[1][i] >= alturaTablero) {
            posicion[1][i] = alturaTablero - 2;
            velocidad[1][i] = - velocidad[1][i];
        }
         
        
    }
    
    free(ax);
    free(ay);
    pthread_exit((void*)0);
}

double factorial(double x) {
    double fac = 1;
    int i;
    for (i=2; i<=x; i++) fac *= i;
    return fac;
}

double mypow(double x, double exp) {
    double result = 1;
    int i;
    
    result = 1;
    for (i=0; i<exp; i++) result *= x;
    return result;
}

double mysin(double alpha) {
    double result = 0;
    int i;
    
    for (i=0; i < MAX_ACCURACY;i++) {
        result += (double)(mypow(-1, i) / factorial((double)(2*i+1))) * (double)mypow(alpha, 2*i+1);
    }
    return result;
}

double mycos(double alpha) {
    double result = 0;
    int i;
    
    for (i=0; i < MAX_ACCURACY;i++) {
        result += (double)(mypow(-1, i) / factorial((double)(2*i))) * (double)mypow(alpha, 2*i);
    }
    return result;
}

